// This code converts LHC Olympics .lhco files into Les Houches Events .lhe files
// written by Roberto Franceschini (2010)
// r2 fix in the definiyion of theta(eta)
// double theta = 2. * atan(exp(- 1. * eta));

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <iterator>
#include <vector>
#include <iostream>
#include <sstream>
#include <string>
#include <cstdlib>
#include "math.h"
#include <stdio.h>

using namespace std;

float ptiso(float HADoverEM) {

	float res;
	res = floor(HADoverEM);
	return res;
}

float etratio(float HADoverEM) {

	float res;
	res = HADoverEM - floor(HADoverEM);
	return res;
}

float HoverEorpTiso(float HADoverEM, int pid) {
	float res;

	if (pid == 13 || pid == -13) {
		res = ptiso(HADoverEM);
	} else {
		res = HADoverEM;
	}

	return res;

}

int LHCOtoPDG(string lhcopid, double ntrack, double btag) {
	int pdgpid = 0;
	if (lhcopid == "0") {
		// photons
		pdgpid = 22;
	}

	if (lhcopid == "1") {
		if (ntrack > 0) {
			//electrons
			pdgpid = -11;
		}

		if (ntrack < 0) {
			// positrons
			pdgpid = 11;
		}
	}

	if (lhcopid == "2") {
		if (ntrack > 0) {
			pdgpid = -13;
		}

		if (ntrack < 0) {
			pdgpid = 13;
		}
	}

	if (lhcopid == "3") {
		if (ntrack > 0) {
			pdgpid = -15;
		}

		if (ntrack < 0) {
			pdgpid = 15;
		}
	}

	if (lhcopid == "4") { //jets
		pdgpid = 211234;

		if (btag == 2) { //tight btag
			pdgpid = 55;
		} else if (btag == 1) { //loose btag
			pdgpid = -55;
		}

	}
	if (lhcopid == "6") { //mET
		pdgpid = 121416;
	}
	return pdgpid;

}

int main(int argc, char* argv[]) {

	if ((argc != 3) and (argc != 2)) {
		cout << " Usage: lhco2lhe input_file output_file" << endl;
		cout << " input_file - input file in LHCO format," << endl;
		cout << " output_file - output file will be in LHE format." << endl;
		return 1;
	}
	cout << " ####################################################### " << endl;

	cout
			<< " This LHCO2LHE revision 2, a code that converts LHCO files into LHE files using particle ID as follows (PDG compliant as much as possible)"
			<< endl;
	cout << " reconstructed photons PID = 22" << endl;
	cout << " reconstructed electron/positron PID = 11/-11" << endl;
	cout << " reconstructed muon/anti-muon PID = 13/-13" << endl;
	cout << " reconstructed tau/anti-tau PID = 15/-15" << endl;
	cout << " reconstructed mET PID = 121416" << endl;
	cout << " reconstructed light flavor jet PID = 211234" << endl;
	cout << " reconstructed b-jet (TIGHT TAG) PID = 55" << endl;
	cout << " reconstructed b-jet (LOOSE TAG) PID = -55" << endl;
	cout << " ####################################################### " << endl;

	int ninout = 0;
	vector < string > lines;
	vector < string > fields;

	string infilename = argv[1];
	string outfilename;

	if (argc == 3) {
		outfilename = argv[2];
		cout << "The given filename for the output is " << outfilename << endl;

	} else if (argc == 2) {
		outfilename = infilename + ".lhe";
		cout << "The autmatic filename for the output is " << outfilename
				<< endl;
	}

	ifstream unw(argv[1]);
	ofstream lhe;

	lhe.open(outfilename.c_str());


	lhe
			<< "<LesHouchesEvents version=\"1.0\">\n<!--\n-->\n<init>\n  2212  2212  0.70000000000E+04  0.70000000000E+04 0 0 10042 10042 3 1\n  0.11750367120E-02  0.10499153982E-04  0.23730090000E-05 661\n</init>"
			<< endl;
	cout
			<< "File Initialized with the following string. Check that this is what you want!"
			<< endl;
	cout
			<< "<LesHouchesEvents version=\"1.0\">\n<!--\n-->\n<init>\n  2212  2212  0.70000000000E+04  0.70000000000E+04 0 0 10042 10042 3 1\n  0.11750367120E-02  0.10499153982E-04  0.23730090000E-05 661\n</init>"
			<< endl;

	cout << "Events stored with LHE informations:" << endl;
	cout << "	0  1.0  0.9118800E+02  0.7816531E-02  0.1300000E+00" << endl;

	lhe.close();
	stringstream eventparticles;
	string temp;
	string tempfield;
	istringstream infield;
	istringstream instream;
	int blockposition = 0;
	string in_pid;
	double eta, phi, pt, jmas, ntrk, btag, HADoverEM;
	int pid;
	double pX, pY, pL, m, E;
	bool write_event = false;

	//LOOP on LINES
	ninout = 0;
	while (getline(unw, temp)) {

		lines.push_back(temp);
		instream.clear();
		instream.str(temp);
		fields.clear();
		while (getline(instream, tempfield, ' ')) {

			if (tempfield != "" && tempfield != " ") {
				fields.push_back(tempfield);
				//cout<<tempfield<<endl;
			}
		}

		// now I know the content of the line
		// cout << "read fields: " << fields.size() << endl;
		//cout << "first field: " << fields[0] << endl;

		if (fields[0].substr(0, 1) == "#") {
			//cout << "was a comment" << endl;
			// this is a comment line, do nothing

		} else if (fields[0] == "0") {
			//this is the trigger line of the event, reset the counter of the particles and
			// ask to write the particles that are waiting to be written
			//cout<<" trigger line"<<endl;

			write_event = true;

		} else {
			//cout<<"Particle line"<<endl;

			// this line is a particle, store it in the container
			ninout++;
			eta = atof(fields[2].c_str());
			phi = atof(fields[3].c_str());
			pt = atof(fields[4].c_str());
			m = atof(fields[5].c_str());
			ntrk = floor(atof(fields[6].c_str())); //some objects have non integer number of tracks
			btag = atof(fields[7].c_str());
			HADoverEM = atof(fields[8].c_str());

			pid = LHCOtoPDG(fields[1], ntrk, btag);

			HoverEorpTiso(HADoverEM, pid);

			double p;
// old and wrong!	double theta = 2 * atan(exp(-2 * eta));
            double theta = 2. * atan(exp(- 1. * eta));

			p = pt / sin(theta);

			pX = pt * sin(phi);
			pY = pt * cos(phi);
			pL = p * cos(theta);

			E = sqrt(pX * pX + pY * pY + pL * pL + m * m);
			if (false) {
				eventparticles << "    " << pid << " 1 " << btag << " " << ntrk
						<< " 0  0 " << pX << " " << pY << " " << pL << " " << E
						<< "  " << m << " " << HoverEorpTiso(HADoverEM, pid)
						<< "  " << etratio(HADoverEM) << endl;
			}
			eventparticles << "    " << pid << " 1 " << btag << " " << ntrk
					<< " 0  0 " << pX << " " << pY << " " << pL << " " << E
					<< "  " << m << " " << "0"/*HoverEorpTiso(
			 HADoverEM, pid)*/<< "  " << "0"/*etratio(HADoverEM)*/<< endl;
			//cout << eventparticles.str()  << endl;
		}

		if ((eventparticles.str()).size() >= 1 && write_event == true) {
			// write the event and reset the flag
			//cout<< "writing event"<<endl;
			lhe.open(outfilename.c_str(), ios::app);
			lhe << "<event>\n" << ninout
					<< " 0  1.0  0.9118800E+02  0.7816531E-02  0.1300000E+00"
					<< endl;
			ninout = 0;
			lhe << eventparticles.str();
			eventparticles.str("");
			//cout<<"now the particle vector is "<<eventparticles.str().size()<<endl;

			lhe << "</event>" << endl;
			lhe.close();

		}

		write_event = false;

	}
	lhe.open(outfilename.c_str(), ios::app);
	lhe << "</LesHouchesEvents>\n" << endl;
	lhe.close();
	return 0;
}

